<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>unsteady-state Flux Balance Analysis (uFBA)</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S4 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S10 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S11 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span style=' font-weight: bold;'>unsteady-state Flux Balance Analysis (uFBA)</span></h1><h2  class = 'S1'><span style=' font-weight: bold;'>James T. Yurkovich</span></h2><div  class = 'S2'><span>Department of Bioengineering and the Bioinformatics and Systems Biology Program, University of California, San Diego USA</span></div><div  class = 'S2'><span>Reviewed by Aarash Bordbar</span></div><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>In this tutorial, we will use unsteady-state Flux Balance Analysis (uFBA) [1] to integrate exo- and endo-metabolomics data [2] into a constraint-based metabolic model for the human red blood cell [3]. The uFBA method allows for bypassing the steady-state assumption for intracellular metabolites that are measured. </span></div><div  class = 'S2'><span>We can model the flux through a metabolic network using a set of linear equations defined by</span></div><div  class = 'S3'><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;block&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;S&lt;/mi&gt;&lt;mo&gt;&amp;sdot;&lt;/mo&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;v&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;b&lt;/mi&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="57" height="18" /></span></div><div  class = 'S2'><span>where </span><span style=' font-weight: bold;'>S</span><span> is the stoichiometric matrix, </span><span style=' font-weight: bold;'>v</span><span> is a vector of fluxes through the chemical reactions defined in </span><span style=' font-weight: bold;'>S</span><span>, and </span><span style=' font-weight: bold;'>b</span><span> represents constraints on the change of metabolite concentrations; at steady-state, </span><span style=' font-weight: bold;'>b </span><span>= 0. If the metabolomics measurements are non-linear (i.e., Fig. 1), then the first step of the uFBA workflow is to identify discrete time intervals which represent linearized metabolic states (Fig. 1). Once discrete states are identified (the raw data if linear), we proceed to estimating metabolite concentration rates of change. For each metabolic state, we can use linear regression to calculate the rate of change of each metabolite concentration. If the rate of change is significant, the model is updated by changing the steady-state constraint from 0 to </span></div><div  class = 'S3'><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;block&quot;&gt;&lt;mrow&gt;&lt;mtable columnalign=&quot;left&quot;&gt;&lt;mtr&gt;&lt;mtd&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;S&lt;/mi&gt;&lt;mo&gt;&amp;sdot;&lt;/mo&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;v&lt;/mi&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;mo&gt;&amp;geq;&lt;/mo&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;b&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/mtd&gt;&lt;/mtr&gt;&lt;mtr&gt;&lt;mtd&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;S&lt;/mi&gt;&lt;mo&gt;&amp;sdot;&lt;/mo&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;v&lt;/mi&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;mo&gt;&amp;leq;&lt;/mo&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;bold&quot;&gt;b&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/mtd&gt;&lt;/mtr&gt;&lt;/mtable&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-17px"><img src="" width="65" height="45" /></span></div><div  class = 'S2'><span>where [</span><span texencoding="\mathbf{b}_1" style="vertical-align:-6px"><img src="" width="16.5" height="20" /></span><span>, </span><span texencoding="\mathbf{b}_2" style="vertical-align:-6px"><img src="" width="16.5" height="20" /></span><span>] represents the 95% confidence interval for each significantly changing metabolite. All unmeasured metabolites are assumed to be at steady-state (i.e., </span><span texencoding="\mathbf{b}_1 = \mathbf{b}_2 = 0" style="vertical-align:-6px"><img src="" width="74.5" height="20" /></span><span>).</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "640.4198730468751" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>                                                    Fig. 1 | Overview of the uFBA workflow.</span></div><h2  class = 'S1'><span>MATERIALS</span></h2><h2  class = 'S1'><span>Equipment Setup</span></h2><div  class = 'S2'><span>Running uFBA requires the installation of a mixed-integer linear programming solver. We have used Gurobi 7.0.0 (http://www.gurobi.com/downloads/download-center) which is freely available for academic use (this workflow has only been tested with Gurobi solvers; use other solvers at your own risk). This tutorial uses the Statistics Toolbox to perform linear regression (if the Statistics Toolbox is not installed, compute linear regression manually; see </span><span style=' font-family: monospace;'>testUFBA.m</span><span>).</span></div><h2  class = 'S4'><span>PROCEDURE </span></h2><h2  class = 'S1'><span>Initialize</span></h2><div  class = 'S2'><span>Running uFBA requires the use of several functions from the COBRA Toolbox.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% no toolbox update, just init</span></span></div></div></div><div  class = 'S2'><span>We first load in sample data. This data is absolutely quantified and has already been volume adjusted such that intracellular and extracellular metabolite concentrations have compatible units. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >tutorialPath = fileparts(which(</span><span style="color: rgb(170, 4, 249);">'tutorial_uFBA.mlx'</span><span >));</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >load([tutorialPath filesep </span><span style="color: rgb(170, 4, 249);">'sample_data.mat'</span><span >]);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% We load the model by readCbModel to make sure it fits to the specifications.</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = readCbModel([tutorialPath filesep </span><span style="color: rgb(170, 4, 249);">'sample_data.mat'</span><span >],</span><span style="color: rgb(170, 4, 249);">'modelName'</span><span >,</span><span style="color: rgb(170, 4, 249);">'model'</span><span >) </span></span></div></div></div><div  class = 'S9'><span>The </span><span style=' font-family: monospace;'>sample_data.mat</span><span> file contains the following variabels:</span></div><ul  class = 'S10'><li  class = 'S11'><span style=' font-family: monospace;'>met_data</span><span>: a matrix containing the exo- and endo-metabolomics data</span></li><li  class = 'S11'><span style=' font-family: monospace;'>met_IDs</span><span>: a cell array containing the BiGG ID for each of the metabolites in </span><span style=' font-family: monospace;'>met_data</span></li><li  class = 'S11'><span style=' font-family: monospace;'>model</span><span>: a modified version [3] of the iAB-RBC-283 COBRA model structure</span></li><li  class = 'S11'><span style=' font-family: monospace;'>time</span><span>: a vector of the time points (in days) at which the metabolite concentrations were measured</span></li><li  class = 'S11'><span style=' font-family: monospace;'>uFBAvariables</span><span>: a struct containing the variables necessary for input into the uFBA algorithm</span></li></ul><div  class = 'S2'><span>In this tutorial, the use of Gurobi is mandatory.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >solverLPOk = changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'gurobi'</span><span >, </span><span style="color: rgb(170, 4, 249);">'LP'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >solverMILPOk = changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'gurobi'</span><span >, </span><span style="color: rgb(170, 4, 249);">'MILP'</span><span >);</span></span></div></div></div><div  class = 'S9'><span></span></div><h2  class = 'S4'><span>Estimate Metabolite Rates of Change (&lt;1 sec.)</span></h2><div  class = 'S2'><span>Next, we run linear regression to find the rate of change for each metabolite concentration.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >changeSlopes = zeros(length(met_IDs), 1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >changeIntervals = zeros(length(met_IDs), 1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(met_IDs)</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >    [tmp1, tmp2] = regress(met_data(:, i), [time ones(length(time), 1)], 0.05);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >    changeSlopes(i, 1) = tmp1(1);</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >    changeIntervals(i, 1) = abs(changeSlopes(i, 1) - tmp2(1));</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S9'><span>The variables </span><span style=' font-family: monospace;'>changeSlopes</span><span> and </span><span style=' font-family: monospace;'>changeIntervals</span><span> contain the metabolite rates of change and 95% confidence intervals, respectively. We will create a new vector, </span><span style=' font-family: monospace;'>ignoreSlopes</span><span>, which contains a 0 for the metabolites whose slopes change significantly and a 1 otherwise:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >tmp1 = changeSlopes - changeIntervals;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >tmp2 = changeSlopes + changeIntervals;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >ignoreSlopes = double(tmp1 &lt; 0 &amp; tmp2 &gt; 0);</span></span></div></div></div><h2  class = 'S4'><span>Integration of Metabolomics Data (&lt;10 min.)</span></h2><div  class = 'S2'><span>Finally, we need to input the data into the uFBA algorithm which is encapsulated in the function </span><span style=' font-family: monospace;'>buildUFBAmodel</span><span>. This function takes as input a COBRA model structure and a struct containing the required input variables (see Table 1). </span></div><div  class = 'S2'><span>Ideally, all metabolites in the model would be measured, resulting in a feasible model. However, experimental limitations limit the number of metabolites that can measured. Thus, when the metabolite constraints are added, the model will most likely not simulate. The uFBA algorithm reconciles the measured metabolomics data and the network structure by parsimoniously allowing unmeasured metabolites concentrations to deviate from steady-state (i.e., </span><span texencoding="\mathbf{S \cdot v = b}" style="vertical-align:-5px"><img src="" width="57.5" height="18" /></span><span>) in order to build a computable model. We refer to the method for deviating unmeasured metabolites from steady-state as "metabolite node relaxation." As part of this procedure, free exchange of extracellular metabolites out of the system is only allowed if (1) the metabolite concentration is measured to be increasing or (2) if the relaxation of a particular extracellular metabolite is required for model feasibility. </span></div><div  class = 'S2'><span>There are five different techniques built into the uFBA method to perform the node relaxation. The technique used in this tutorial is an MILP optimization that minimizes the number of unmeasured metabolites relaxaed from steady-state; this choice effectively minimizes the changes made to the model in order to achieve feasibilitiy. Full details for this and all other node relaxation techniques can be found in [1]. Sinks are added for each of the relaxed metabolite nodes, and the flux through each of these sinks is minimized while still allowing the model to simulate. The minimimum value is then multiplied by a relaxation factor </span><span style=' font-family: monospace;'>lambda</span><span> (Table 1) and used as the bound for the sink reaction.</span></div><div  class = 'S2'><span>Full details for the algorithm are provided in the original publication [1]. </span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "664.716796875" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>                                       Table 1 | Inputs and outputs of the </span><span style=' font-family: monospace;'>buildUFBAmodel</span><span> function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >uFBAvariables.metNames = met_IDs;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >uFBAvariables.changeSlopes = changeSlopes;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >uFBAvariables.changeIntervals = changeIntervals;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >uFBAvariables.ignoreSlopes = ignoreSlopes;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >uFBAoutput = buildUFBAmodel(model, uFBAvariables);</span></span></div></div></div><div  class = 'S9'><span>The output contains the resulting model (</span><span style=' font-family: monospace;'>uFBAoutput.model</span><span>):</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >model_ufba = optimizeCbModel(uFBAoutput.model)</span></span></div></div></div><h2  class = 'S4'><span>References</span></h2><div  class = 'S2'><span>[1] A Bordbar*, JT Yurkovich*, G Paglia, O Rolfsson, O Sigurjonsson, and BO Palsson. "Elucidating dynamic metabolic physiology through network integration of quantitative time-course metabolomics." </span><span style=' font-style: italic;'>Sci. Rep. </span><span>(2017). doi:10.1038/srep46249. (* denotes equal contribution)</span></div><div  class = 'S2'><span>[2] A Bordbar, PI Johansson, G Paglia, SJ Harrison, K Wichuk, M Magnusdottir, S Valgeirsdottir, M Gybel-Brask, SR Ostrowski, S Palsson, O Rolfsson, OE Sigurjonsson, MB Hansen, S Gudmundsson, and BO Palsson. "Identified metabolic signature for assessing red blood cell unit quality is associated with endothelial damage markers and clinical outcomes." </span><span style=' font-style: italic;'>Transfusion</span><span> (2016). doi:10.1111/trf.13460.</span></div><div  class = 'S2'><span>[3] A Bordbar, D McCloskey, DC Zielinski, N Sonnenschein, N Jamshidi, and BO Palsson. "Personalized Whole-Cell Kinetic Models of Metabolism for Discovery in Genomics and Pharmacodynamics." </span><span style=' font-style: italic;'>Cell Systems</span><span> (2015). doi:10.1016/j.cels.2015.10.003.</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% *unsteady-state Flux Balance Analysis (uFBA)*
%% *James T. Yurkovich*
% Department of Bioengineering and the Bioinformatics and Systems Biology Program, 
% University of California, San Diego USA
% 
% Reviewed by Aarash Bordbar
%% INTRODUCTION
% In this tutorial, we will use unsteady-state Flux Balance Analysis (uFBA) 
% [1] to integrate exo- and endo-metabolomics data [2] into a constraint-based 
% metabolic model for the human red blood cell [3]. The uFBA method allows for 
% bypassing the steady-state assumption for intracellular metabolites that are 
% measured. 
% 
% We can model the flux through a metabolic network using a set of linear equations 
% defined by
% 
% $$\mathbf{S}\cdot \mathbf{v}=\mathbf{b}$$
% 
% where *S* is the stoichiometric matrix, *v* is a vector of fluxes through 
% the chemical reactions defined in *S*, and *b* represents constraints on the 
% change of metabolite concentrations; at steady-state, *b* = 0. If the metabolomics 
% measurements are non-linear (i.e., Fig. 1), then the first step of the uFBA 
% workflow is to identify discrete time intervals which represent linearized metabolic 
% states (Fig. 1). Once discrete states are identified (the raw data if linear), 
% we proceed to estimating metabolite concentration rates of change. For each 
% metabolic state, we can use linear regression to calculate the rate of change 
% of each metabolite concentration. If the rate of change is significant, the 
% model is updated by changing the steady-state constraint from 0 to 
% 
% $$\begin{array}{l}\mathbf{S}\cdot \mathbf{v}\;\ge {\mathbf{b}}_1 \\\mathbf{S}\cdot 
% \mathbf{v}\;\le {\mathbf{b}}_2 \end{array}$$
% 
% where [$\mathbf{b}_1$, $\mathbf{b}_2$] represents the 95% confidence interval 
% for each significantly changing metabolite. All unmeasured metabolites are assumed 
% to be at steady-state (i.e., $\mathbf{b}_1 = \mathbf{b}_2 = 0$).
% 
% 
% 
% Fig. 1 | Overview of the uFBA workflow.
%% MATERIALS
%% Equipment Setup
% Running uFBA requires the installation of a mixed-integer linear programming 
% solver. We have used Gurobi 7.0.0 (http://www.gurobi.com/downloads/download-center) 
% which is freely available for academic use (this workflow has only been tested 
% with Gurobi solvers; use other solvers at your own risk). This tutorial uses 
% the Statistics Toolbox to perform linear regression (if the Statistics Toolbox 
% is not installed, compute linear regression manually; see |testUFBA.m|).
%% PROCEDURE 
%% Initialize
% Running uFBA requires the use of several functions from the COBRA Toolbox.

initCobraToolbox(false) % no toolbox update, just init
%% 
% We first load in sample data. This data is absolutely quantified and has already 
% been volume adjusted such that intracellular and extracellular metabolite concentrations 
% have compatible units. 

tutorialPath = fileparts(which('tutorial_uFBA.mlx'));
load([tutorialPath filesep 'sample_data.mat']);
% We load the model by readCbModel to make sure it fits to the specifications.
model = readCbModel([tutorialPath filesep 'sample_data.mat'],'modelName','model') 
%% 
% The |sample_data.mat| file contains the following variabels:
%% 
% * |met_data|: a matrix containing the exo- and endo-metabolomics data
% * |met_IDs|: a cell array containing the BiGG ID for each of the metabolites 
% in |met_data|
% * |model|: a modified version [3] of the iAB-RBC-283 COBRA model structure
% * |time|: a vector of the time points (in days) at which the metabolite concentrations 
% were measured
% * |uFBAvariables|: a struct containing the variables necessary for input into 
% the uFBA algorithm
%% 
% In this tutorial, the use of Gurobi is mandatory.

solverLPOk = changeCobraSolver('gurobi', 'LP');
solverMILPOk = changeCobraSolver('gurobi', 'MILP');
%% 
% 
%% Estimate Metabolite Rates of Change (<1 sec.)
% Next, we run linear regression to find the rate of change for each metabolite 
% concentration.

changeSlopes = zeros(length(met_IDs), 1);
changeIntervals = zeros(length(met_IDs), 1);
for i = 1:length(met_IDs)
    [tmp1, tmp2] = regress(met_data(:, i), [time ones(length(time), 1)], 0.05);
    changeSlopes(i, 1) = tmp1(1);
    changeIntervals(i, 1) = abs(changeSlopes(i, 1) - tmp2(1));
end
%% 
% The variables |changeSlopes| and |changeIntervals| contain the metabolite 
% rates of change and 95% confidence intervals, respectively. We will create a 
% new vector, |ignoreSlopes|, which contains a 0 for the metabolites whose slopes 
% change significantly and a 1 otherwise:

tmp1 = changeSlopes - changeIntervals;
tmp2 = changeSlopes + changeIntervals;
ignoreSlopes = double(tmp1 < 0 & tmp2 > 0);
%% Integration of Metabolomics Data (<10 min.)
% Finally, we need to input the data into the uFBA algorithm which is encapsulated 
% in the function |buildUFBAmodel|. This function takes as input a COBRA model 
% structure and a struct containing the required input variables (see Table 1). 
% 
% Ideally, all metabolites in the model would be measured, resulting in a feasible 
% model. However, experimental limitations limit the number of metabolites that 
% can measured. Thus, when the metabolite constraints are added, the model will 
% most likely not simulate. The uFBA algorithm reconciles the measured metabolomics 
% data and the network structure by parsimoniously allowing unmeasured metabolites 
% concentrations to deviate from steady-state (i.e., $\mathbf{S \cdot v = b}$) 
% in order to build a computable model. We refer to the method for deviating unmeasured 
% metabolites from steady-state as "metabolite node relaxation." As part of this 
% procedure, free exchange of extracellular metabolites out of the system is only 
% allowed if (1) the metabolite concentration is measured to be increasing or 
% (2) if the relaxation of a particular extracellular metabolite is required for 
% model feasibility. 
% 
% There are five different techniques built into the uFBA method to perform 
% the node relaxation. The technique used in this tutorial is an MILP optimization 
% that minimizes the number of unmeasured metabolites relaxaed from steady-state; 
% this choice effectively minimizes the changes made to the model in order to 
% achieve feasibilitiy. Full details for this and all other node relaxation techniques 
% can be found in [1]. Sinks are added for each of the relaxed metabolite nodes, 
% and the flux through each of these sinks is minimized while still allowing the 
% model to simulate. The minimimum value is then multiplied by a relaxation factor 
% |lambda| (Table 1) and used as the bound for the sink reaction.
% 
% Full details for the algorithm are provided in the original publication [1]. 
% 
% 
% 
% Table 1 | Inputs and outputs of the |buildUFBAmodel| function.

uFBAvariables.metNames = met_IDs;
uFBAvariables.changeSlopes = changeSlopes;
uFBAvariables.changeIntervals = changeIntervals;
uFBAvariables.ignoreSlopes = ignoreSlopes;

uFBAoutput = buildUFBAmodel(model, uFBAvariables);
%% 
% The output contains the resulting model (|uFBAoutput.model|):

model_ufba = optimizeCbModel(uFBAoutput.model)
%% References
% [1] A Bordbar*, JT Yurkovich*, G Paglia, O Rolfsson, O Sigurjonsson, and BO 
% Palsson. "Elucidating dynamic metabolic physiology through network integration 
% of quantitative time-course metabolomics." _Sci. Rep._ (2017). doi:10.1038/srep46249. 
% (* denotes equal contribution)
% 
% [2] A Bordbar, PI Johansson, G Paglia, SJ Harrison, K Wichuk, M Magnusdottir, 
% S Valgeirsdottir, M Gybel-Brask, SR Ostrowski, S Palsson, O Rolfsson, OE Sigurjonsson, 
% MB Hansen, S Gudmundsson, and BO Palsson. "Identified metabolic signature for 
% assessing red blood cell unit quality is associated with endothelial damage 
% markers and clinical outcomes." _Transfusion_ (2016). doi:10.1111/trf.13460.
% 
% [3] A Bordbar, D McCloskey, DC Zielinski, N Sonnenschein, N Jamshidi, and 
% BO Palsson. "Personalized Whole-Cell Kinetic Models of Metabolism for Discovery 
% in Genomics and Pharmacodynamics." _Cell Systems_ (2015). doi:10.1016/j.cels.2015.10.003.
##### SOURCE END #####
-->
</div></body></html>